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Abstract 

We present a phenomenological theory of the interplay between nematic order and supercon- 
ductivity in the vicinity of a vortex induced by an apphed magnetic field. Nematic order can be 
strongly enhanced in the vortex core. As a result, the vortex cores become elliptical in shape. For 
the case where there is weak bulk nematic order at zero magnetic field, the field-induced eccentricity 
of the vortex core has a slow power-law decay away from the core. Conversely, if the nematic order 
is field-induced, then the eccentricity is confined to the vortex core. We discuss the relevance of 
our results to recent scanning tunneling microscopy experiments on FeSe (Song et al, Science 332, 
1410 (2011)). 
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I. INTRODUCTION 



The unconventional superconductors have a rich phase diagram determined by the in- 
terplay of multiple competing, or coexisting, types of order. Nematic order (which breaks 
the C4 symmetry of the underlying square lattice down to C2) has been shown to emerge 
in certain regimes of the phase diagrams of the copper-oxide [1-6] and the iron-based [7-13] 
superconductors. In the latter case, the nematic order accompanies (and in some cases, 
precedes) the magnetic order which occurs at a wavevector that breaks the lattice rotational 
symmetry. 

Recently, the structure of the vortex cores in the mixed state of clean FeSe films was 
studied by means of scanning tunneling microscopy (STM) [14]. Strong anisotropy was 
observed in the zero bias conductance map around the cores, which have an eccentricity of the 
order of unity. Although the lattice structure of FeSe at low temperature is orthorhombic[15], 
it has been claimed [14] that the crystalline anisotropy (of the order of a few tenths of a 
percent) is too small to explain the large anisotropy of the vortex cores, which is likely to 
have an electronic origin. 

This experiment raises several questions, some of which we address in this paper: assum- 
ing that there is an electronic nematic order in superconducting FeSe, what is its microscopic 
origin? What is its relation to superconductivity - e.g., are these two types of order com- 
peting? Is the nematic order localized in the vortex cores (and hence stabilized by the 
application of the magnetic field), or does it extend throughout the system (and is only 
apparent in the STM spectrum near the cores)? 

Here, we study the structure of the vortex core using a phenomenological Landau- 
Ginzburg (LG) theory in terms of two competing order parameters. Using our LG analysis 
we have calculated the structure of an isolated vortex in the presence of the nematic order. 
Our main result is that by looking at the profile of the gap near the vortex core, it is pos- 
sible to distinguish between two different configurations of the nematic order, namely the 
presence of a localized nematic order within the superconducting vortex as opposed to the 
presence of a long range nematic order in the system. If the nematic order is localized at the 
core, the superconducting gap should be anisotropic only near the core and the anisotropy 
decays exponentially as we move away from the core. On the other hand, if the nematic 
order is long-ranged, the superconducting gap should exhibit an anisotropy which decays as 
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a power law. If the nematic order is near its critical point, there is a large region in which the 
anisotropy of the gap depends logarithmically on the distance, eventually crossing over to a 
power law. Moreover, we find qualitative differences in the shape of the contours of constant 
gap around the core in the different cases. If the nematic order exists only in the cores, 
the equal-gap contours tend to be elliptical; if the nematic order is long-ranged, we find 
that the gap function tends to develop a "four-lobe" structure, with more pronounced higher 
harmonics. These features can be sought in STM experiments by mapping the magnitude 
of the gap around the core as a function of position. 

The paper is organized as follows: In section II we introduce the LG functional with the 
two competing order parameters and carry out a preliminary analysis in the absence of the 
anisotropy. In section III, we investigate the mean-field phase diagram of a single vortex. 
In section IV, we introduce the anisotropy and perform a numerical minimization of the 
functional, commenting on the interesting features. Finally, in section V, we present our 
analytical results explaining the various interesting features observed by minimizing the free 
energy. 

II. MODEL 

We consider a LG type free energy for two competing order parameters: a complex field 
describing the superconducting order parameter, and a real field </>, which describes a 
nematic order that competes with the superconducting order parameter. The form of the 
free energy density is given by 
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Apart from the standard free energy contributions arising due to and we have a com- 
petition term, controlled by 7 (> 0), and a term that gives rise to different effective masses 
for \l/ in the two directions, which is controlled by Ai. J-" is invariant under a rotation by 90 
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degrees, represented by 



X ^ y, 
y -X, 

^ -0. (5) 

We will be interested in the limit of A — oo, where A is the London penetration depth, so 
that we can neglect the coupling to the electromagnetic field. At the outset, we set A2 = 0, 
since the A2 term is small compared to the Ai term in the limit where is small. It is 
convenient to define the coherence length of and the healing length of (p as 

Taking the unit of distance to be l^, we can recast the above free energy in a more transparent 
form as follows. 
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where / = /^//^,7s = 7V'o/0O'-^ = ^i/'^^l'^o^ = x/l^,y/l^, ip = '^/ipo, and = 
From now on, we will drop the tilde symbols. 

For A 7^ 0, a short-distance cutoff has to be imposed on Eq. 7. Otherwise, the system is 
unstable towards developing modulations of ip with sufficiently short wavelength. We discuss 
the instability in Appendix A. In practice, we will mostly ignore this issue, assuming that 
there is a short-distance cutoff (which is provided by the finite grid used in our numerical 
calculations) . 

Before we begin our analysis, let us comment about the choice of parametrization in this 
problem. We would like to think of this problem in terms of a fixed 7 < 1. Then on choosing 
a particular ratio of the length scales of and ip, we still have one degree of freedom left in 
terms of the masses or the stiffnesses of the two order parameters, which is fixed by tuning 

Is- 

If we assume that 7^ > 1, then the uniform ground state is given by -0 = 1 and = 0. 
This also constrains to be localized around the vortex cores, by making the mass term for 



positive deep inside the superconducting region. If we further assume that the nematic 
order is small, such that <^ tp^, then we can essentially ignore the feedback of (p on ijj. 
Therefore, we will first find the full profile of = \E'o, the isolated vortex solution, in the 
absence of the nematic order and use that to find the form of the nematic order. Then \l/o 
satisfies the following asymptotic relations: 



^o(p) 



*o(p) ~ Clpe' 



2 \lp 



p->r^ 
P<^r^ 



(8) 
(9) 



where p = r/l^, r being the radius in the original coordinate system, and C is a dimensionless 
constant. In general, it is difficult to find the solution of the full LG equation for \&o for all 
p analytically. Therefore we obtain the vortex solution \l/o = f{lp)e^^ for all p by minimizing 
the functional in Eqn. 7 numerically in the absence of 0. 

The numerical solution conforms to the two asymptotic expressions above. It is interesting 
to note that \&o does not recover from the vortex core to its bulk value exponentially, but 
rather as a power law [16, 17]. The behavior of 0(p) in the vicinity of a vortex with A = 
was studied by Ref. 16. 



III. PHASE DIAGRAM 



We will now describe the mean-field phase diagram of a single vortex in the presence of 
a competing nematic order. There are three possible phases: in phase I, </> = everywhere; 
in phase II, </> vanishes at large distance from the vortex core but becomes non-zero near 
the vortex core due to the suppression of the competing ip field to zero at the core; and in 
phase III, 7^ even far away from the core. A non-zero solution for is favored whenever 
the smallest eigenvalue e of the following eigenvalue problem [16]: 



V?-l + 7.[/(^P)]^ 



0(x) = e0(a:), (10) 



satisfies e < 0. In order to find the phase diagram, we solve this eigenvalue problem nu- 
merically on a discrete grid. The boundary between phases I and II is the locus of points 
at which the smallest eigenvalue satisfies e = 0. For 7^ < 1, becomes long-ranged, corre- 
sponding to phase III. The resulting phase diagram is shown in Fig. 1. This phase diagram 
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Figure 1: The phase diagram in the (7^, I) plane obtained by solving Eqn. 10 numerically on a grid 
with n = 100 and Ax = 1. The regions with qualitatively different solutions for (f) are marked. 
Phase I has no nematic order with = everywhere, phase II has nematic order localized in the 
vortex core, and phase III has long-range nematic order away from the core. The blue squares 
correspond to points which we explore in more detail later. The dashed line 7^ = 1 represents the 
boundary between phases II and III. 

is strictly valid as long as 7^ > 7. If this is not the case, then the state with uniform nematic 
background and no superconductivity is energetically favorable over any other state. 

The physics behind the phase diagram can be understood as follows. When ^ 
we are forcing the nematic order to coexist with superconductivity in a large region. 
This is unfavorable energetically due to the competition term 70^|\&p. Therefore, when 
7s > 1, there is no 7^ solution. If 7^ < 1, becomes non-zero even far away from 
the vortex core. In the opposite limit of <^ l^, the nematic order exists deep within 
the superconducting vortex. Since there is very little overlap between the two order 
parameters, the system can afford to have a higher value of critical 7^ below which there is a 
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nontrivial nematic order. This explains the increasing trend of the critical 7^ for decreasing /. 

In the <^ case, it is possible to give an analytical expression for the phase boundary 
between regions I and 11. The equation for this curve is given by, 

^ 4(C7F' ^^^^ 

where C is the constant which appears in Eqn. 9. The details of this computation are 
diskussed in Section VA. 

We are now in a position to include the effect of the anisotropy and investigate the 
structure of the vortex cores in different regions of the phase diagram described above. 



IV. VORTEX PROFILE IN THE DIFFERENT REGIMES 

We now turn to diskuss the characteristics of the vortex profile in the different regimes 
shown in Fig. 1. To solve for the vortex profile, we minimize the free energy (7) with 
respect to ip and (p numerically on a disk geometry. This is equivalent to solving the coupled 
Landau-Ginzburg equations with Neumann boundary conditions, as we diskuss in Appendix 
B. Many of the features found in the numerical solution can be understood analytically, as 
we diskuss in the next section. 

We can expand both tp and (p in terms of the different angular momentum channels 
e*"^). The term proportional to A only couples angular momentum channels that differ 
by 2 units of angular momentum in ip. Therefore, in the presence of 0, the bare vortex 
solution (~ e*^) gives rise to components of the form e^*^, e~*^, etc. Similarly, the feedback 
of the superconducting order on (p gives rise to the generation of the even harmonics, i.e. $0 
gives rise to terms proportional to e^*^ and e~^*^. It is also possible to have a solution with 
only the even harmonics oiip, in which case, the vortex is absent. These two solutions do not 
mix with each other and therefore we shall focus on the solution in the presence of the vortex. 

In light of this, we expand the order parameters as 

^P{p, 9) = Y^ ^n(p)e^('"+')', 0(p, 0) = J2 Mpy'''', n G integers (12) 

n n 
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In terms of the expansions in Eqn. 12, the free energy density can be written as, 
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and we are interested in minimizing J pdpFp. We shall minimize the above free energy for 
a given system size and for only a fixed number of harmonics at a time. We have kept n 
harmonics for and 0, where for any given n (odd) we take all the harmonics ^ _i to 
i = {n — l)/2, and similarly for (p. We have tried n = 3, 5 and found no substantial qualitative 
change in the results that we shall quote here, indicating that the results converge even with 
only 3 harmonics. We consider a system on a disk of radius p = 100. 

Below, we describe the results in regions II and III of the phase diagram, and on the 
critical line dividing them (In region I, where there is no nematic order, we get the regular 
circularly symmetric vortex) . The specific values of 7^ , / which were used are marked by 
blue squares in the phase diagram in Fig. 1. 



A. Region II 

In this region, we expect to obtain a solution with a non-zero uniform {tpl away from the 
vortex core and a non-zero (p localized near the vortex core, decaying exponentially away 
from the core (given that 7 < 7s and 7^ > 1). The contour plot for l^/^p is shown in Fig. 2. 
The parameters used here are 7s = 2.0, / = 0.1, 7 = 1.0, A = 20.0. 

As can be seen in the figure, the core has an elliptical shape because of the interaction 
with the nematic order which coexists with superconductivity in the core region. As we go 
away from the core, the contours of equal become more and more isotropic, due to the 
rapid decay of the nematic order away from the core. 
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Figure 2: Nematic order in phase II. The profiles for the different order parameters for 7<j = 2.0, 1 = 
0.1,7 = 1.0, A = 20.0 for a system size of 100 {N = 100, Ax = 1). (a) Harmonics of ip (sohd) and 
$0 (dashed) (b) Contour plot of IV'P- The superconducting coherence length is 10, in units of 

B. Region III 

This region in the phase diagram corresponds to the case where there is a uniform nematic 
background coexisting with superconductivity, even away from the vortex core. In this 
regime, as we move away from the core, (p goes to a constant and tp remains anisotropic. In 
Fig. 3, the harmonics and \l/„i are almost constant for large p. Moreover, \E'i = — \E'„i 
for large p. The contour plot of reveals a large anisotropic "halo" around the core, with 
a non-elliptical shape. 

Far away from the core, where (p is constant, the Landau- Ginzburg equations can be 
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Figure 3: Nematic order in phase III. The profiles for the different order parameters for 7^ = 0.5, 1 = 
0.5,7 = 0.4, A = 0.5 for a system size of 100 {N = 100, Ax = 1). (a) Harmonics of (solid) and 
$0 (dashed) (b) Contour plot of The superconducting coherence length is 2, in units of Ifj^. 

solved analytically, showing that the anisotropy in {ipl'^ decays as a power law in this case. 
We diskuss this solution in the next section. 



C. Critical case 



Finally, we diskuss the critical line separating regions II and III in Fig. 1, in which 
the field is critical far away from the core. Naively, one would expect to go as 1/p 
asymptotically in this regime. However, depending on the details of the solution at small p, 
there may be an intermediate regime in which ~ In p, eventually crossing over to 1/p at a 
larger distance. This feature is diskussed in more detail in the next section. Fig. 4a shows 
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Figure 4: Nematic order at the critical point between phase III and phase I. The profiles for 
the different order parameters for 7^ = 1.0,/ = 3.0,7 = 0.9, A = 0.07 for a system size of 100 
(N = 100, Ax = 1). (a) Harmonics of ip (solid) and <I>o (dashed) (b) Contour plot of IV'P- The 
superconducting coherence length is 1/3, in units of Z^. 

■ip and (j) in the critical regime, with 7^ = 1.0,/ = 3.0,7 = 0-9 ^^^d A = 0.07. Indeed, we 
observe that decays slowly away from the core. "^±1 also have long tails. The contour plot 
for I'lpl'^ shares features that are similar to the behavior in region III, namely a long-range, 
non-elliptical anisotropic halo. It is shown in Fig. 4b. 

In the next section, we analyze the asymptotic behavior of the solution in the critical 
case, showing that the anisotropic component of falls off as ~ (Alnp/p2)cos(2^) at 
intermediate p, crossing over to ~ (A/p^) cos(20) at sufficiently large p. 
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V. ANALYTICAL TREATMENT 



In this section, we propose various analytical arguments to explain the different features 
that were observed above by carrying out the minimization numerically. In subsection VA, 
we diskuss the solution for the nematic order in the presence of superconductivity. Then in 
subsection VB, we analyze region III of the phase diagram (Fig. 1). Finally, in subsection 
VC, we study the linearized GL equations in order to explain some of the other interesting 
features that were observed earlier. 



A. Phases of the nematic order 



Here we will briefly review the solution for (f), and supplement it with some further details. 
The LG-equation for 4>{p), assuming that A = 0, is given by. 



l+ls[f{lp)? + <l>' 



0, 



(14) 



We shall now be interested in solving the linearized version of the above equation, which is 
justified for 7^ > 1. For p <^ this becomes equivalent to solving the problem 
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This is identical to solving the Schrodinger equation for the 2D quantum harmonic oscillator. 
We know that 0'(p = 0) = 0. The solution for this equation is given by 



-^f^pV2 ^aWlsCn^P^ 



(16) 



where Cn{x) are the Laguerre polynomials. The profiles of 0(p) for a few different values of 
7^(7^/^ are shown in Fig. 5. 

At this point, we can also describe how we obtained the equation for the phase boundary 
between regions I and II in the phase diagram (Fig. 1). In this case, (p is non-zero only very 
close to the center of the core. We can therefore expand ip around p = and keep only the 
leading order term (Eq. 9). Eqn.15 can be re- written as, 

i2;2 



2 2 ^ 



e)0. 



(17) 



Non-trivial solutions exist for e < 0. The above equation is the Schrodinger equation for 
a quantum harmonic oscillator in two dimensions with m = l,a;^ = 7s(C/)^. Then the 
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Figure 5: The profiles of (/> as a function of p for different values of 7<jC^/^ over a distance of one 
correlation length of (f), l^j,. 



smallest eigenvalue which corresponds to the zero point energy of the oscillator leads to the 
following equation for the curve 



1 



Is 



A{Clf 
in the limit of small /. 

On the other hand, for p ^ we have to solve 

v^- 1 + 7/1 



(18) 



0, 



(19) 



from which we see that (j){p) ~ p~ ' exp(— ^^7^ — Ip) . However, there is a fine-tuned point 
at 7s = 1, at which the field is critical far away from the vortex core. The full equation 
for (h becomes 



1 



+ 



0, 



(20) 



</'(p) = a/I + l'^ / (Ip)- Note that the 1/p solution can be obtained only for specific boundary 
conditions. For generic boundary conditions, with — )■ as p — 00, the solution is neverthe- 



less asymptotic to \/T+P/{lp) at large p. If 0(1) ^ 1, for instance, then 0(p) ~ Aln(po/p) 
at intermediate values of p, where A and po are constants. 0(p) crosses over to 0(p) ~ 1/p 
at radii of the order of p* ~ 1/A (see Appendix C). This behavior reflects itself in the 
asymptotic decay of the anisotropy of the field ip away from the vortex core, as we saw in 
Sec. IV. 
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B. Coexistence of superconductivity and nematic order 



In this subsection, we are interested in analyzing region III of the phase diagram, in which 
superconductivity and nematicity coexist even far away from the core. Let us assume, for 
simphcity, that far away from the core can be replaced by a constant. The effect of a 
constant is to render the effective masses in the two directions different. Therefore, if we 
re-scale the coordinates as 



X 



y' = (21) 

where a = 2A0/^, then this problem now becomes identical to the isotropic problem we had 
solved in the beginning of section II. The solution for \l'o can then be written in terms of 
the new coordinates as, 

f(r') 



m = c(l --—,], c=Jl-^ (22) 



Note that due to the presence of the background nematic order, \&o does not tend to 1 
asymptotically. We now go back to our original coordinate system x, y by expanding the 
above result to linear order in a. Then we get. 

In the above expression, the first bracket corresponds to \E'o) the second bracket corresponds 
to while the last one represents \E'i. It is interesting to observe that asymptotically, 
and — approach the same constant value. We observe this feature in Fig. (3a). However, 
the harmonics do not recover to their asymptotic value as a power law, which is a result 
of the boundary conditions that were imposed while minimizing the free energy in the disk 
geometry (see Appendix B). 

From Eqn. 23, we can evaluate the form of l^/^p and find that, 

= c^(l - ^) - ^cos(2^) + 0(a^) (24) 
Therefore, asymptotically, lip]"^ is isotropic and the anisotropy decays as ~ a cos(26')/r^. 

14 



C. Linearized GL analysis 



In this section, we shall carry out an analysis of the linearized LG equations, to give an 
analytical explanation for some of the features that we have observed by carrying out the 
full minimization. For the sake of simplicity, let us ignore the feedback on resulting in 
the generation of the higher harmonics and assume that is isotropic (i.e. 0(p, 6') = 0(p)). 
Then, the linearized LG equations for the harmonics of il) can be written as. 



-A0(p) 



dl + {An + b)^ + ^ '-\ '-]^n+i{p) 



p p 

dp bn-l(p) + 9, + — )vl/„+i(p) 



(25) 



P 

There are some features of the problem that cannot be deduced from a study of the linearized 
version of the problem, which include the overall scale and sign of and the signs of the 
different harmonics of ip. 

In the limit of p ^ i.e. inside the vortex core, at leading order \l/n(p) ~ p"", where 
a = |2n+ 1|. This is a necessary condition for the harmonics to be well behaved in the limit 
of p -> 0. 

On the other hand, in the limit of p ^ /~^, i.e deep inside the superconducting region, the 
homogenous solution for the above equation gives exponentially damped solutions for all the 
\l/„^0; i-6. the anisotropy is short ranged. Moreover, the source term, which is proportional 
to A0 and is itself exponentially damped (Region II), is also not strong enough to give rise 
to any long ranged solution. 

However, when is critical (i.e. 7^ = 1), the source term leads to the presence of long 
tails in the harmonics. In the regime where falls off logarithmically while \l/o is a constant, 
at leading order just follow 0, i.e. they also fall off logarithmically (with prefactors of 
equal magnitude but opposite sign) and have a correction of the form In p/p^. On the other 
hand, when crosses over to the power law form, at leading order "$±1 also fall off as ±l/p 
with a correction of order 1/p'^. 
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VI. CONCLUSION 



We have studied the interplay between nematic order and superconductivity in the pres- 
ence of a vortex. If the nematic order coexists with superconductivity in the vicinity of a 
vortex core, the couphng between the two order parameters leads to an elongated shape 
of the core. We diskuss two distinct scenarios: in one the nematic order coexists with su- 
perconductivity everywhere (i.e., even far away from the vortex core), whereas in the other 
the competition between the two order parameter suppresses the nematic order in the bulk, 
and nematicity only exists close to the core where the superconducting order parameter is 
diminished. Both scenarios lead to an anisotropic core. However, we show that they can, 
in principle, be distinguished by the way the anisotropy of the superconducting gap decays 
away from the core. If the nematicity exists only near the core, the anisotropy in the su- 
perconducting gap decays exponentially; if it exists throughout the sample, we expect the 
gap anisotropy to decay as where r is the distance from the core. Moreover, there are 

qualitative differences in the shape of the core in the two cases. In the former case, in which 
only the core region is nematic, the contours of equal gap tend to be more or less elliptical. 
In the latter case, the contours of equal gap tend to develop non-elliptical shapes with a 
four-petal pattern. Therefore, analyzing the gap profiles measured by STM around a vortex 
could reveal the nature of the nematic ordering - whether it is localized at the vortex core, 
or coexists with superconductivity in the bulk. 

So far, we have diskussed the structure of an isolated vortex at the mean-field level. 
However, if the nematic ordering is favored only within a vortex core, an isolated vortex 
cannot have static nematic order, since either thermal or quantum fluctuations would destroy 
such order. Static nematic order is only possible when the density of vortices is finite. The 
coupling between the nematic halos of different vortices scales as J^s ~ exp[— d/ — '^s^4>)], 
where d ~ 1/\/B is the inter- vortex distance (B is the applied magnetic field). The system 
can be described by an effective two-dimensional transverse field Ising model with a spin-spin 
interaction Jcr and a 5-independent transverse field. (Note that, unlike Ref. [16], we are 
considering a thin film, rather than a three-dimensional system.) This model has a nematic 
transition at a certain critical B, which should be seen, e.g., by measuring the anisotropy of 
the vortex cores as a function of B. If an external rotational symmetry breaking field exists, 
as is presumably the case in FeSe due to the small orthorhombic lattice distortion[15], the 
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electronic nematic transition is smoothed out. However, one still expects a sharp crossover 
as a function of magnetic field if the orthorhombic distortion is sufficiently weak. 

The microscopic origin of the anisotropic vortex cores observed in FeSe[l4] remains to 
be understood. It is likely that it originates from electronic nematicity rather than from 
the lattice distortion, since the experimentally reported orthorhombic distortion seems too 
small to produce such a large effect. The electronic nematic order could have an orbital 
character[ll, 18-20]. Alternatively, it could arise from a field-induced magnetic ordering[21] 
at a wavevector (tt, 0) or (0, tt) in the one iron unit cell, which is necessarily accompanied by a 
nematic component (similar to the ordering in the iron arsenides). Although static ordering 
of this type has not been observed in the iron selenides[22], it remains to be seen if they 
develop a static ordering in the presence of an applied magnetic field. Neutron scattering 
experiments revealed a magnetic resonance at this wavevector in the superconducting state 
of FeTeSe[23]. Moreover, ordering at such wavevector s nearly nests the electron and hole 
pockets, and therefore it is expected to couple strongly to superconductivity, explaining why 
the resulting anisotropy of the vortex cores is so large. 

Note added: After this work was submitted for publication, another manuscript [24] that 
studied the experimental features observed in FeSe [14] came to our attention. In this 
paper, the authors study the effect of orbital ordering on the vortex structure in a two band 
model, by solving the Bogoliubov-de Gennes equations. This study is complementary to our 
phenomenological Ginzburg-Landau approach. 
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Appendix A: Instabilities of the free energy 

An interesting feature associated with the LG functional introduced in section II is that 
there is an instability to a state with modulated This arises due to a competition between 
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two terms in the free energy, namely the (j) and ifP' terms. Let us suppose that does not 
vary spatially and if) = Pe^'^^. Then at leading order, the contribution to the free energy 
from (J) is of the form 

-^^= f^-7?^V' + V/3V (Al) 



V 2% 272^ 

From the above expression, we see that for a sufficiently large Ag^, it becomes energetically 
favorable to gain energy from the second term by condensing a large negative value of (j). 

By extremizing the above with respect to 0, we obtain 0^ = — A/3^g^/ ( — -h] ■ Hence, 

the contribution to free energy from 0m is oc — /3^A^g^. This energy gain from a non-zero 
q always dominates over the energy cost of order for a sufficiently large q. In order to 
prevent this instability, we have to add a term of the form (\V'^ip\^ /2l'^ to the free energy, 
which is an allowed term from the underlying symmetry of the problem. We now want to 
obtain some restrictions on (. 

First of all, ( should be such that it prevents the instability. This gives us a lower bound on 
the value of (. At the same time, ( should be small enough so that it should not change the 
physics significantly. This gives us an upper bound on the value of Therefore, we obtain, 

AS^/^ 

^ < C « 1 (A2) 

The above expression is not valid when becomes critical, i.e. when 7^ = 1, /3 = 1. In this 
case, we have to compare with 0^. 

However, when we minimized the free energy in section IV, we did not have to include the 
above term with a finite ( as for a sufficiently small A, the cutoff in q arising from the discrete 
lattice prevented this instability from showing up. 



Appendix B: Effect of boundary terms 

In general, when we derive the GL equations from the free energy, there is a surface term 
arising from the gradient terms in the energy which can be ignored in the limit of an infinite 
system size. However, for a finite sized system, the boundary term does play an important 
role. Let us consider only the contribution of the gradient term of the superconducting order 
parameter in the free energy, in the absence of any nematic order. Then we have, 

Fgrad = J rfV|V^|^ (Bl) 
Ptot -^grad ~l~ -^local) 

(B2) 
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where Fiocai contains the usual l^/^p, terms. On varying ip* by Sif)* in Ftoti we obtain for 
a finite system (up to other variations due to Fiocai denoted by ...), 

- f d^rSip*V^ip + f Sij*{Vip) ■nds + ... = 0, (B3) 

J J surface 

where fids is the area element, normal to the boundary. When we solve for ip in the interior of 
the region, only the first term contributes and the boundary term can be ignored. However, 
when we solve for ip on the boundary, only the surface term plays a role, since it can be 
thought of as appearing with an infinite weight of the form J dr6{r — R) where R is the 
radius of the disk on which we are minimizing the free energy and is the Dirac-delta 

function. Therefore, in order to solve for tp, we have to solve for —V'^ip + ... = in the 
interior of the region subject to the boundary condition Vip ■fi\r=R = (Neumann boundary 
conditions) . 

Now in the presence of a constant nematic background {(po), the gradient term in the free 
energy is, 

Fgrad = J dxdy {l + a)\d^ip\^ + il-a)\dytp\'' , (B4) 

where a = 2X(pQp. As we did earlier, on carrying out the variation over tp* this amounts to 
solving for —{1 + a)d'^ip — {1 — a)dyip + ... = subject to the boundary condition, Dip-fi = 0, 
where 

Dip = (^{l + a)d^^,{l-a)dy^^, fi = {cos9,sm9). (B5) 
In polar coordinates, this condition can be written as. 



drip + a 



cos 29 drip — — — ^OqiP 
r 



0. (B6) 



These boundary conditions mean, in particular, that the current perpendicular to the 
boundary is zero. In our numerical calculations, we have used a disk geometry; therefore the 
boundaries are found to have a significant effect whenever we are considering a non-circularly 
symmetric solution, in particular in the regime where the nematic order is non-zero even 
far away from the core. We circumvent this problem, however, by taking a sufficiently large 
system and considering the solution only close to the vortex core, where the boundary effects 
are small. 
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Figure 6: Numerical solutions of Eq. CI. a) Solutions with boundary conditions (j)(1000) = and 
with various values of 0(1), on a log-log scale, b) Solution with boundary condition = 0.125, 
0(1000) = 0, on a semilog scale. The dashed line is a fit to the form A + B log{p) for small p. 

Appendix C: Asymptotics of in the critical case 

In this appendix, we analyze the asymptotics of the field far away from the vortex core 
in the case 7^ = 1, in which the nematic order is critical. In this case, and for p ^ 1, the 
Landau-Ginzburg equation for (Eq. 14) becomes 



0. 



(CI) 
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This non-linear equation admits the solution (j){p) = 1/p [16]. This solution is valid, however, 
for specific initial conditions, e.g., 0(1) = 1, = ~1- Physically, the initial conditions for 
Eq. CI are determined by the details of the vortex profile at short distances, determined by 
Eq. 14. Nevertheless, one can make some general statements about the asymptotic behavior 
of the solution. If, for some arbitrary po such that po ^ 1/^ (far from the core), (p satisfies 
(po) ^ 1/PO) then it is justified to neglect the 0^ term in Eq. CI. Then, the solution 
close to Po behaves as (p) ^ A — B In (p/po), where A, B are determined by the initial 
conditions. This can only be valid, however, up to a point at which (p*) ~ 1/p*, i-e., 
at distances which are much smaller than the length scale set by the initial condition of Eq. 
CI. At longer distances, we expect a crossover to (p) ^ 1/p. 

In Fig. 6, we present a numerical solution of Eq. CI with boundary conditions 0(1000) = 
and various values for 0(1). When 0(1) = 1, we get 0(p) ~ 1/p (where the deviations 
are due to the boundary condition at p = 1000). For smaller 0(1), there is an intermediate 
region where does not follow a power law, eventually crossing over to 1/p at larger p. 0(p) 
is approximately logarithmic in the intermediate region, as shown in Fig. 6b. 

Physically, we expect that < 1 (since 0=1 corresponds to the equilibrium value of 
in the absence of superconductivity). Therefore, there is an intermediate logarithmic region, 
which becomes parametrically large in the limit of small 0. 
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